function output = DetermineOptimalL4andL5State(IC,T,u,DU,VU,TU,dt)
options=odeset('RelTol',1e-13,'AbsTol',1e-13);
options_fmincon = optimoptions('fmincon', 'MaxFunctionEvaluations', 20000, ...
    'MaxIterations', 20000,'ConstraintTolerance',1e-6,"Algorithm","interior-point",...
    "EnableFeasibilityMode",true,...
    "SubproblemAlgorithm","cg");
%

to = 0:dt:T;
Xo = IC;
[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S'; t = t';

S_i = zeros(size(S));
for ii = 1:length(S)
    S_i(:,ii) = C2I_primary(S(:,ii),u,DU,VU,t(ii));
end

ind = find(S_i(3,:)==max(S_i(3,:)));

S_z_max = S_i(:,ind);

InclinationRotationAngle = 7.25;
RAANRotationAngle = 73.67+0.014*(2023-1850);
PoleDirection = Rotation([0;0;1],RAANRotationAngle,'Degrees')*...
                Rotation([1;0;0],InclinationRotationAngle,'Degrees')*...
                [0;0;1];

S_z_max = [S_z_max(1);S_z_max(2);0];
PoleDirect = [PoleDirection(1);PoleDirection(2);0];
alpha = acosd(dot(S_z_max,PoleDirect)/(norm(S_z_max)*norm(PoleDirect)));

T_mat = Rotation([0;0;1],alpha,'Degrees');
r_i = zeros(3,length(S));
for ii = 1:length(S)
    r_i(:,ii) = T_mat*S_i(1:3,ii);
end

% Initial condition is the decending node.

%we know that the Earth's position is 60 degress behind L4.

alpha = acos(dot(S_z_max,PoleDirect)/(norm(S_z_max)*norm(PoleDirect)));

to = 0:-0.0001:(-alpha);
Xo = IC;
[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S'; t = t';

to = 0:dt:T;
Xo = S(:,end);
[t,S_r_L4] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S_r_L4 = S_r_L4'; t = t';
S_i_L4 = zeros(size(S_r_L4));
for ii = 1:length(S_r_L4)
    S_i_L4(:,ii) = C2I_primary(S_r_L4(:,ii),u,DU,VU,t(ii));
end

epoch = datetime(2023,1,1,0,0,0) + seconds(t*TU);

output.S_i = S_i_L4;
output.S_r = S_r_L4;
output.epoch = epoch;

%%

% % OpenCRTBP_u(u); hold on
% plot3(S_r_L4(1,:),S_r_L4(2,:),S_r_L4(3,:),'k-','LineWidth',2);
% xlabel('AU');
% ylabel('AU');
% zlabel('AU');
% 
% ind_min = find(S_r_L4(3,:)==min(S_r_L4(3,:)));
% ind_max = find(S_r_L4(3,:)==max(S_r_L4(3,:)));
% 
% plot3(S_r_L4(1,ind_min),S_r_L4(2,ind_min),S_r_L4(3,ind_min),'rd','MarkerSize',8);
% timeepochtext = epoch(ind_min);
% text(S_r_L4(1,ind_min),S_r_L4(2,ind_min),S_r_L4(3,ind_min), datestr(timeepochtext,'mm/dd'),'FontSize',15,'BackgroundColor','none','Color','r')
% 
% plot3(S_r_L4(1,ind_max),S_r_L4(2,ind_max),S_r_L4(3,ind_max),'rd','MarkerSize',8);
% timeepochtext = epoch(ind_max);
% text(S_r_L4(1,ind_max),S_r_L4(2,ind_max),S_r_L4(3,ind_max), datestr(timeepochtext,'mm/dd'),'FontSize',15,'BackgroundColor','none','Color','r')
% 
% % Find index near half of up and down.
% 
% zamphalf = S_r_L4(3,ind_max)/2;
% temp_S_r_z = find(abs(S_r_L4(3,:)-zamphalf)==min(abs(S_r_L4(3,:)-zamphalf)));
% plot3(S_r_L4(1,temp_S_r_z),S_r_L4(2,temp_S_r_z),S_r_L4(3,temp_S_r_z),'rd','MarkerSize',8);
% timeepochtext = epoch(temp_S_r_z);
% text(S_r_L4(1,temp_S_r_z),S_r_L4(2,temp_S_r_z),S_r_L4(3,temp_S_r_z), datestr(timeepochtext,'mm/dd'),'FontSize',15,'BackgroundColor','none','Color','r')
% 
% zamphalf = S_r_L4(3,ind_min)/2;
% temp_S_r_z = find(abs(S_r_L4(3,:)-zamphalf)==min(abs(S_r_L4(3,:)-zamphalf)));
% plot3(S_r_L4(1,temp_S_r_z),S_r_L4(2,temp_S_r_z),S_r_L4(3,temp_S_r_z),'rd','MarkerSize',8);
% timeepochtext = epoch(temp_S_r_z);
% text(S_r_L4(1,temp_S_r_z),S_r_L4(2,temp_S_r_z),S_r_L4(3,temp_S_r_z), datestr(timeepochtext,'mm/dd'),'FontSize',15,'BackgroundColor','none','Color','r')


% xlim([0.44,0.54])
% ylim([0.8,0.9])
% zlim([-0.3,0.3])

% L5
end




























